log using log_main.smcl, replace

********************************************************************************
*
*	The Political Geography of the January 6 Insurrectionists
*   PS: Political Science & Politics
*   Replication
*
*   Robert A. Pape
*   Kyle D. Larson
*   Keven G. Ruby
*
*	Create all tables and figures in article
*   
*   v2024-01-08
*
********************************************************************************

* User Installed Packages

/*

Jann, Ben. 2023. "ESTOUT: Stata Module to Make Regression Tables." Statistical 
Software Components, February. https://ideas.repec.org//c/boc/bocode/s439301.html.

Pisati, Maurizio. 2018. "SPMAP: Stata Module to Visualize Spatial Data." 
Statistical Software Components, January. https://ideas.repec.org//c/boc/bocode/s456812.html.

Winter, Nick. 2021. "COMBOMARGINSPLOT: Stata Module to Combine the Saved Results
from Multiple Calls to Margins into One Marginsplot." Statistical Software 
Components, November. https://ideas.repec.org//c/boc/bocode/s457804.html.

*/

/*
* spmap
ssc install spmap
* combomarginsplot
ssc install combomarginsplot
* esttab
ssc install estout
*/

* load the data
use PS_Jan6countydata.dta, clear


* standardize explanatory variables and controls to mean = 0 and sd = 1

foreach var of varlist wpdecline mfgempdecline trump2020vromney pctnhwhite2020 urban distdc countypop {
	
	egen z`var' = std(`var')
	
} 

////////////////////////////////////////////////////////////////////////////////
// Table 1. Counties that Fought for Trump vs. Counties that Voted for Trump  //

* Column 1: Counties with >= 1 Insurrectionists

table () insurrectbinary, ///
stat(mean urban mfgempdecline wpdecline pctnhwhite2020 pop_tot_2020 distdc) ///
nformat(%9.2g)

// Column 2: Counties where Trump 2020 > Romney 2012	

gen tvrbinary = cond(trump2020vromney>0,1,0)
 
table () tvrbinary, ///
stat(mean urban mfgempdecline wpdecline pctnhwhite2020 pop_tot_2020 distdc) ///
nformat(%9.2g)

	


//////////////////////////////////////////////////////////////////////////////// 
// Figure 1. Geographic Distribution of Indicted Insurrectionists             //

spmap insurrectcount  using "us_county_coor_proj.dta", id(_ID) ///
clmethod(custom) clbreaks(0, 0.9, 2, 5, 10, 20) ///
fcolor(black*.05 red*.20 red*.5 red*1.3 red*2) osize(.1pt ..) ocolor(gs10 ..) ///
legorder(lohi)  ///
legend(pos(3) ring(1) col(1) rowgap(0.3) size(*1.5) margin(l = 0.6 b = -.01) ///
    symysize(*0.9)  alignment(top)  ///
    label(2 "0 (2,683)")    ///
    label(3 "1 to 2 (355)")    ///
    label(4 "3 to 5 (80)")   ///
    label(5 "6 to 10 (18)")   ///
    label(6 "> 10 (5)")   ///
    title("{bf:County Count of}" "{bf:Insurrectionists}", pos(11) size(small)) ///
	note("Mean All: 0.3""Mean > 0: 2.1", margin(t=2))) ///
	polygon(data("us_state_coor_proj.dta") fcolor(none)  ///
        ocolor(white%70 ..) osize(vthin ..)) 

	graph export figure1.tif, width(1000) replace



////////////////////////////////////////////////////////////////////////////////
// Table 2. Modeling Insurrection Propensity of Counties                      //
  
 
 * Negative binomial regressions
 	 
	// Model 1: White Population Decline 
	nbreg insurrectcount 							///
	zwpdecline 										///
	zpctnhwhite2020 urban zdistdc 					///
	zcountypop	 									///
	, robust irr

	estimates store m1

	// Model 2: Manufacturing Loss
	nbreg insurrectcount 							///
	zmfgempdecline 									///
	zpctnhwhite2020 urban zdistdc					///
	zcountypop	 									///
	, robust irr

	estimates store m2

	// Model 3: Populist Support for Trump
	nbreg insurrectcount 							///
	ztrump2020vromney 								///
	zpctnhwhite2020 urban zdistdc 					///
	zcountypop	 									///
	, robust irr

	estimates store m3

	// Model 4: Full
	nbreg insurrectcount 							///
	zwpdecline zmfgempdecline ztrump2020vromney		///
	zpctnhwhite2020 urban zdistdc					///
	zcountypop	 									///
	, robust irr	

	estimates store m4

	// Model 5: Full, Counties w/ Pop >0 and <1 SD
	nbreg insurrectcount 							///
	zwpdecline zmfgempdecline ztrump2020vromney		///
	zpctnhwhite2020 urban zdistdc					///
	zcountypop	 									///
	if zcountypop > 0 & zcountypop < 1 				///
	, robust irr	

	estimates store m5

	// Model 6: Full, Counties > Median Pop
	nbreg insurrectcount 							///
	zwpdecline zmfgempdecline ztrump2020vromney		///
	zpctnhwhite2020 urban zdistdc					///
	zcountypop	 									///
	if pop_tot_2020 > 25750 						/// 
	, robust irr	

	estimates store m6
	
	// Model 7: Full, Counties w/ Pop <3 SD
	nbreg insurrectcount 							///
	zwpdecline zmfgempdecline ztrump2020vromney		///
	zpctnhwhite2020 urban zdistdc					///
	zcountypop	 									///
	if zcountypop < 3 								/// 
	, robust irr	

	estimates store m7	


* Export models 1-7 as HTML table

esttab m1 m2 m3 m4 m5 m6 m7 using Table2.html 					///
, replace se compress aic bic eform constant 					///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) 	///
title("Table 2. Modeling Insurrection Propensity of Counties" ) ///
drop(lnalpha) 													///
coeflabels(														///		 
		zwpdecline			"White Population Decline, 2010-2020" ///
		zmfgempdecline		"Decline in Manufacturing Employment, 1970-2020" ///
		ztrump2020vromney	"Populist Support for Trump"	///			
		zpctnhwhite2020		"% Population Non-Hispanic white" ///
		urban				"Urban = 1"							///
		zdistdc				"Distance to DC (1,000 kms)" ///
		zcountypop			"County Population (100,000s)" ///
		_cons 				"Constant"	/// 
			) 	///
order(zwpdecline zmfgempdecline ztrump2020vromney zpctnhwhite2020 urban zdistdc ///
zcountypop _cons) lines nogaps wide nodepvar ///
note("Note: Negative binomial regression. Reporting IRRs. All continuous covariates standardized." ///
"Robust standard errors in parentheses.")


////////////////////////////////////////////////////////////////////////////////
// Figure 2. Effect on predicted county count of insurrectionists             //

* caclulate margins for models 1-7 
foreach num of numlist 1 4 5 6 7 {
	
	di "`num'"
	estimates restore m`num'
	margins, at(zwpdecline=(-3(1)3)) atmeans saving(z`num'w, replace)	
}

// (a) White Population Decline
combomarginsplot  z1w z4w,  noci  legend(ring(1) pos(6) colfirst col(1)) ///
graphregion(margin(vsmall)) scheme(plotplain) ///
addplot(hist zwpdecline if zwpdecline >= -3 & zwpdecline <= 3, width(1) ///
legend(order(1 "Model 1" 2 "Model 4")) ///
title("{bf:(a) White Population}" "{bf:Decline}") ///
ytitle("Predicted Count of Insurrectionists", size(medium)) ///
xtitle(" ", size(small)) ///
yaxis(2) yscale(alt axis(2)) percent   ///
ylabel(0(.1).5, labcolor(black*.9) axis(1)) ///
xlabel(-3(1)3) ///
xscale(range(-3 3)) ///
xlabel(-3 "-3" -2 "-2" -1 "-1" 0 "0" 1 "1" 2 "2" 3 "3", labcolor(black*.9)) ///
ylabel(, labcolor(black*.9) axis(2))  ///
ytitle(" ", axis(2) orientation(rvertical))  fcolor(white%10) lcolor(black*.5)) ///
xsize(6.5) ysize(4.5) name(z1, replace)
graph export z1.png, width(1000) replace


// (b) Manufacturing Loss
foreach num of numlist 2 4 5 6 7 {
	
	di "`num'"
	estimates restore m`num'
	margins, at(zmfgempdecline=(-3(1)3)) atmeans saving(z`num'm, replace)	
}

combomarginsplot  z2m z4m,  noci  legend(ring(1) pos(6) colfirst col(1)) ///
graphregion(margin(vsmall)) scheme(plotplain) ///
addplot(hist zmfgempdecline if zmfgempdecline >= -3 & zmfgempdecline <=3, width(1) ///
legend(order(1 "Model 2" 2 "Model 4")) ///
title("{bf:(b) Manufacturing}""{bf:Decline}") ytitle(" ", size(medium)) ///
xtitle("Standard Deviations", size(small)) ///
yaxis(2) yscale(alt axis(2)) percent   ///
ylabel(0(.1).5, labcolor(black*.9) axis(1)) ///
xlabel(-3(1)3) ///
xlabel(-3 "-3" -2 "-2" -1 "-1" 0 "0" 1 "1" 2 "2" 3 "3", labcolor(black*.9)) ///
ylabel(, labcolor(black*.9) axis(2))  ///
ytitle(" ", axis(2) orientation(rvertical))  fcolor(white%10) lcolor(black*.5)) ///
xsize(6.5) ysize(4.5) name(z2, replace)
graph export z2.png, width(1000) replace


// (c) Populist Support for Trump
foreach num of numlist 3 4 5 6 7 {
	
	di "`num'"
	estimates restore m`num'
	margins, at(ztrump2020vromney=(-3(1)3)) atmeans saving(z`num't, replace)	
}

combomarginsplot  z3t z4t,  noci  legend(ring(1) pos(6) colfirst col(1)) ///
graphregion(margin(vsmall)) scheme(plotplain) ///
addplot(hist ztrump2020vromney if ztrump2020vromney>=-3 & ztrump2020vromney<=3, width(1) ///
legend(order(1 "Model 3" 2 "Model 4")) title("{bf:(c) Populist Support}""{bf:for Trump}") ytitle(" ", size(medium)) ///
xtitle(" ", size(small)) ///
yaxis(2) yscale(alt axis(2)) percent   ///
ylabel(0(.1).5, labcolor(black*.9) axis(1)) ///
xlabel(-3(1)3) ///
xlabel(-3 "-3" -2 "-2" -1 "-1" 0 "0" 1 "1" 2 "2" 3 "3", labcolor(black*.9)) ///
ylabel(, labcolor(black*.9) axis(2))  ///
ytitle("Percent of Sample", axis(2) orientation(rvertical))  fcolor(white%10) lcolor(black*.5)) ///
xsize(6.5) ysize(4.5) name(z3, replace )
graph export z3.png, width(1000) replace


graph combine z1 z2 z3, ycommon rows(1) scheme(plotplain)
graph export figure2.tif, replace width(3000)


* Clean-up intermediate files

foreach num of numlist 4/7 {
	
	erase z`num'w.dta
	erase z`num'm.dta
	erase z`num't.dta
}

foreach num of numlist 1/3 {
	
	erase z`num'.png
}

erase z1w.dta
erase z2m.dta
erase z3t.dta



////////////////////////////////////////////////////////////////////////////////
// Table 3. Effect of Political Circumstances and White Population Decline    //

* Negative binomial regressions

	// Model with Competitive Interaction 
	nbreg insurrectcount 			/// 
	c.zwpdecline##i.competitive 	///
	zpctnhwhite2020 urban zdistdc 	///
	zcountypop	 					///
	, robust irr

	estimates store m8

	// Model with Outbidding Interaction 
	nbreg insurrectcount 			/// 
	c.zwpdecline##i.objected 		///
	zpctnhwhite2020 urban zdistdc							///
	zcountypop	 					///
	, robust irr

	estimates store m9
	
	// Model with Binary Measure of Frustration 
	nbreg insurrectcount 			/// 
	c.zwpdecline##i.BidenWonBy20 	///
	zpctnhwhite2020 urban zdistdc 							///
	zcountypop	 					///
	, robust irr

	estimates store m10

	
* Export models 8-10 as HTML table	
	
esttab m8 m9 m10 using Table3.html 								///
, replace se compress aic bic eform constant 					///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) 	///
title("Table 3. Effect of Political Circumstances and White Population Decline") ///
drop(_cons lnalpha zpctnhwhite2020 urban zdistdc zcountypop ) ///
coeflabels(	///		 
		_cons "Constant"														/// 
		zwpdecline 					"White Population Decline, 2010-2020" ///
		1.competitive 				"Contested County"	///
		1.objected					"House Rep Objected = 1" ///
		1.BidenWonBy20				"Biden Won by 20%" 	///
		1.competitive#c.zwpdecline 	"Contested County * White Pop Decline" ///	
		1.objected#c.zwpdecline 	"House Rep Objected * White Pop Decline" ///
		1.BidenWonBy20#c.zwpdecline "Biden Won by 20% * White Pop Decline" 	///	
		) 	///
order(zwpdecline 1.competitive 1.competitive#c.zwpdecline 1.objected 1.objected#c.zwpdecline ///
1.BidenWonBy20 1.BidenWonBy20#c.zwpdecline ) 			///
lines nogaps wide 	///
nodepvars 			///
note("Note: Negative binomial regression. Reporting IRRs. All continuous covariates standardized." ///
"Robust standard errors in parentheses.")



log close 
translate  log_main.smcl log_main.pdf 
